Evidence for the role of fluctuations in the thermodynamics of nanoscale drops and 
the implications in computations of the surface tension 
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Test area deformations are used to analyse vapour-liquid interfaces of Lennard- Jones particles 
by molecular dynamics simulation. For planar vapour-liquid interfaces the change in free energy is 
captured by the average of the corresponding change in energy, the leading-order contribution. This 
is consistent with the commonly used mechanical (pressure tensor) route for the surface tension. 
By contrast for liquid drops one finds a large second-order contribution associated with fluctuations 
in energy. Both the first- and second-order terms make comparable contributions, invalidating the 
mechanical relation for the surface tension of small drops. The latter is seen to increase above the 
planar value for drop radii of ~ 8 particle diameters, followed by an apparent weak maximum and 
slow decay to the planar limit, consistent with a small negative Tolman length. 



It is striking that though almost a century has passed 
since Gibbs formulated his thermodynamic theory of 
curved interfaces, there is still widespread controversy 
about the dependence of the surface tension on the cur- 
vature (size of a drop), and the validity of the mechan- 
ical route to the surface tension-^ 3 . The formal ap- 
proach of Gibbs is intimately connected with the rela- 
tions of Laplace, Ap = 2j s /R s , and Tolman, 7(i?)/7oo = 
1 — 2(5 oc /i?+- ■ ■ , for drops of radius R. Here, Ap = pi—p g 
is the pressure difference inside (1) and outside (g) the 
drop, 7 S = l(Rs) is the interfacial tension associated 
with the surface of tension R s , the value for the pla- 
nar gas-liquid surface, and the Tolman^ length Soo is de- 
fined relative to the radius of the equimolar surface R e 
as = hm fls ^ 00 (i? e - R s ). 

There are 3 basic routes to the definition of the ten- 
sioni: thermodynamic (Gibbs and Tolman), mechanical 
(Laplace and Young) , and statistical mechanical (density 
functional and related theories). The thermodynamic 
and mechanical routes are macroscopic theories so there 
has been much debate about their applicability to small 
systems such as nanoscale liquid drops or bubbles. One 
key question is whether the mechanical relations based on 
the pressure virial (formulated in terms of the appropri- 
ate tensorial components) that make use of the concept of 
the bulk pressure of the coexisting states are appropriate 
at these length scales for curved surfaces. 

While the Laplace equation essentially defines the ratio 
7 s /i? s , the first-order form of Tolman's theory is appro- 
priate only for sufficiently large drops. One can view 
8oo as the leading-order correction to the tension of a 
planar surface. Despite its fundamental role in studies 
of interfacial properties of curved surfaces and theories 
of nucleation, there is still much controversy as to even 
the sign of Soo- Microscopic statistical mechanical ap- 
proaches including square gradient theories (SGTs)^^, 
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curvature expansions of the planar interface^, and den- 
sity functional theories (DFTs), including local^— and 
non-locaii^— treatments, have led to conflicting views 
on the magnitude and sign of S^, as well as the cur- 
vature dependence of the surface tension. The widely 
accepted view from this body of work is that Soo < 
and that there is a small maximum in js(R) as the drop 
radius is decreased, then followed by a sharp decrease. 
This is supported by studies on the penetrable sphere 
modell^ (which can be solved exactly at the mean-field 
level at zero temperature) where one finds a negative Tol- 
man length (Soo = ~~ (7 /2), with a the molecular diameter. 

By contrast, the vast majority of computer simulation 
studies suggest that <5oo > 0. In most simulations of liq- 
uid drops, the mechanical route to the interfacial tension 
is employed, usually involving an integration of the gradi- 
ent of the normal component of the pressure tensor from 
the centre of the drop to the bulk vapour phased—. In 
this case one predicts a monotonous decrease in the sur- 
face tension with increasing curvature (decreasing drop 
radius) from the planar limit (infinite radius); this would 
correspond to Soo > throughout. As was pointed out 
early on by Schofield and Henderson 2 - there are funda- 
mental problems in employing local pressure tensors and 
the associated definition of the internal pressure for mi- 
croscopic (high curvature) drops. This leads to a mis- 
match in the free energy of the formation of a drop deter- 
mined via the mechanical and thermodynamic routes as 
observed in simulation^. Macroscopic thermodynamic 
routes based on a combination of the Laplace and Tolman 
relations have been employed^ but also suffer from the 
ill-definition of the internal pressure and density of the 
liquid. One can also estimate the interfacial tension from 
the free energy change accompanying a volume deforma- 
tion of spherical surfaces using a virial-like expression 21 ; 
these results for the surface tension are in disagreement 
with those obtained from the direct mechanical route. 
Recent grand canonical simulations^, and a thermody- 
namic analysis of large drops based on Laplace- Tolman 
relations 2 ^ both now appear to suggest a small negative 
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5oo, which is consistent with the findings of DFT. 

The aim of this Letter is to use a new method for the 
calculation of the surface tension of small liquid drops in 
molecular simulation, highlighting the role played by the 
fluctuations in the energy of deformation. The method 
relies on the thermodynamic definition of the surface ten- 
sion and is thus free from the inconsistencies associated 
with the application of the mechanical route. A variant of 
the test-area (TA) method 2 ^ is used where small virtual 
perturbations are made in the box dimensions of systems 
with interfaces to obtain the change in free energy associ- 
ated with the corresponding change in surface area. For 
a fluid drop of radius R the change in the Hclmholtz free 
energy F is expressed thcrmodynamically as 1 
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with S the entropy, V g j the vapour and liquid volumes, T 
the temperature, \x the chemical potential, N the number 
of particles, A the interface area, and C the conjugate 
variable for R. The surface tension of a drop is given by 
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where the minimal interfacial tension 7 S defines R s and 
corresponds to taking C = 0. The change in free energy 
AF due to a virtual change in area AA can be expressed 
as the average of the Boltzmann factor of the correspond- 
ing change in configurational energy AU— 



AF 



- fc Tm(exp(-f )) (3) 
= (AU) - ^{(AC/ 2 ) - (AC/) 2 } (4) 
+ ^{< A ^V3<At/ 2 )(At/)+2(At/) 3 }. 

The averages are over configurations of the unperturbed 
reference system. In Eq. (0| AF is expressed as 
a perturbation series to 0((At/ 3 )), where the first- 
order average of the change in energy is AFi = (AU), 
the second-order energy fluctuation term is AF2 — 

-{(At/ 2 ) - (A[/) 2 }/(2fcT), and the third-order contri- 
bution is denoted by AF3. The full Boltzmann form, 
Eq. is employed in, e.g., the test- particle approach 
for the chemical potential 25 , or the volume perturbation 
method for the pressure 2 ^ and the pressure tensor—. 

The tension is obtained as the change in free energy per 
unit area for infinitesimal perturbations to O ((At/ 3 )): 
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The leading term, AFi = (At/), corresponds to the me- 
chanical work involved in changing the area of the inter- 
face, which can be directly associated with the so-called 
virial expression for the tension 2 ^ (expressed in terms of 
averages of the appropriate components a of the virial, 
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FIG. 1: Test-area deformations of a planar liquid-vapour in- 
terface of the LJ-TS fluid. MD simulations of N = 749 par- 
ticles in a periodic box of d 1II16I1S10I1S Ijx — L y = 7.885a and 
Lz = 6L X at T* = kT/e = 0.8 over 3 x 10 6 timesteps. The de- 
formations correspond to changes in the box dimensions (par- 
ticle coordinates) of L' x = L x y/1 + AA* , L' y = L y ^l + AA* , 
and L' z = L a /(1 + AA*). a) The contributions A Fi /AA and 
AF2 /AA to the change in free energy per unit area (in units of 
e/a 2 ) (AA* < 0, +; AA* > 0, X; average, •). The interfacial 
tension 7* = 7<r 2 /e is obtained by extrapolation to AA* = 0. 
b) The distribution V(AU) of the change in energy (relative 
to its average in units of e) scaled at the maximum peak height 
for different relative deformations AA* . The width (standard 
deviation, ctau) is depicted in the inset. 



(x a (dll/dxa)), at the Hookian linear-response level). 
The corresponding entropic contribution due to the de- 
formation is 28 : TAS = {(U) (AU) - (UAU)} /{kT). 

In the case of the a planar interface, it is well known 
that the interfacial tension can be obtained formally from 
the virial expressions^, i.e., entirely from the leading- 
order contribution of Eq. ([5])- This is exemplified for 
a planar vapour-liquid interface of Lennard- Jones par- 
ticles (of diameter a and well depth e, truncated and 
shifted TS at r c = 2.5a) as shown in Fig. 1. A pla- 
nar interface is first stabilised during an NVT molecular 
dynamics (MD) simulation of the inhomogeneous system 
with a liquid slab in the centre of the box separated by 
two vapour regions. The change in configurational en- 
ergy due to small test changes in the dimensions of the 
box such that the area is increased or decreased at fixed 
overall volume is then computed to estimate the various 
contributions in Eq. ([5]) ; the limit of infinitesimal defor- 
mations is obtained by extrapolation to AA — > 0. From 
Fig. la it is clear that only the leading mechanical term 
in AFi/AA contributes to the interfacial tension of a 
planar interface, confirming the validity of the pressure- 
tensor route in this case. The fluctuation term AF2/AA 
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is very small by comparison and does not contribute to 
the tension in an appreciable way; this is also true for 
the third-order term. In Fig. lb we plot the distribution 
V(AU) of the change in configurational energy (relative 
to (At/)) for different area perturbations; the distribu- 
tion is well represented by a Gaussian, the width of which 
(cau) decreases to zero with A A — s> 0, consistent with a 
very small AF 2 /AA - 1 x lQ- & e/a 2 . 

The overall physical picture is fundamentally different 
for a nanosized spherical drop of liquid in contact with 
its vapour. Once the drop has been stabilised its size can 
be characterised from the density profile p(r) as a func- 
tion of the distance r from its centre by calculating the 
Gibbs dividing surface = (p v — pi)^ 1 J ' dr r 3 dp(r)/dr, 
corresponding to an area of A = AttR 2 . Virtual pertur- 
bations from the equilibrium spherical drop geometry are 
made with test changes in the dimensions of the simula- 
tion cube: two of the Cartesian axes are decreased (or 
increased) in length and the third is increased (or de- 
creased) such that the overall volume remains constant. 
The perturbed states correspond to ellipsoidal drops of 
prolate (or oblate) shape which always have larger sur- 
face areas than the original drop, A A > 0. This essen- 
tially corresponds to the longest P2 (Legendre polyno- 
mial) capillary- wave oscillations possible for the drop 3 ; 
the capillary-wave surface tension is equivalent to the 
thermodynamic one at least to leading order in curvature 
0(1/ R). Averages are then accumulated over very long 
runs of ~ 1.5 x 10 9 timesteps, corresponding to microsec- 
ond runs for typical molecular parameters. The term 
AFi/AA is more than two orders of magnitude larger 
in the case of the drop than for a planar interface sys- 
tems of comparable size (cf. Fig. la and Fig. 2a). The 
most significant difference is the large contribution from 
the second-order energy "fluctuation" term AF2/ AA for 
the drop, which was negligible for the planar interface; 
this term is now comparable in magnitude but of oppo- 
site sign to the first-order term. The third-order terms 
remains essentially negligible. Thus both the first- and 
second-order terms contribute to the surface tension of 
the drop. A thermodynamic characteristic of the drop 
is thus the non-vanishing (and large) fluctuation term, 
which is clearly an indication of an additional entropic 
contribution. This can be seen in the distribution of the 
change in configurational energy for different TA pertur- 
bations (Fig. 2b). The data are again well described as 
Gaussians, but now the width does not vanish in the limit 
of infinitesimal deformations, i.e., liniAA-s-o &AU 7^ 
and correspondingly limAA->o AF 2 ^ 0. The fact that 
liniAA^o AF3 ~ for both the planar and curved sys- 
tems suggests symmetrical Gaussians. 

The dependence of the surface tension computed from 
ellipsoidal deformations as a function of the drop size 
(for systems with N — 749 to 11,334) is depicted in 
Fig. 3. Here the tension is computed for R e rather than 
R s , though 7 e = 7 S to 0(1/R 2 ). The behaviour ob- 
tained with our thermodynamic TA approach does not 
support the findings obtained from a standard pressure- 
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FIG. 2: Test-area ellipsoidal deformations of a spherical drop 
of LJ-TS liquid of radius (Re) = 5.55a in coexistence with 
its vapour. MD simulations of N = 749 particles in a peri- 
odic box of dimensions L x = L y — L z = 20a at T* — 0.8 
over 1.5 x 10 9 timesteps. The deformations correspond to 
changes in the box dimensions (particle coordinates) of 14 = 
LiVl + AA*, L'j = LjvT + AA*, and L' k = L k /(l + AA") 
(where i, j, k denotes any of the Cartesian axes), a) The con- 
tributions AFi/AA, AF 2 /AA, and AF 3 /AA to the change 
in free energy per unit area (in units of e/er 2 ) (prolate, +; 
oblate, X; average, •). The tension 7* = 7C" 2 /e is obtained 
by extrapolation to AA* = 0. b) The distribution V(AU) 
of the change in energy scaled at the maximum peak height 
for different relative deformations AA*; the width gau is de- 
picted in the inset. The Gaussians for the planar interface are 
shown dotted (note the very small scale in comparison). 
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FIG. 3: The surface tension of spherical drops of LJ-STS flu- 
ids with average radii (R e /a) = 5.55, 6.12, 8.04, 11.7 and 
14.6 at T* = 0.8 from TA ellipsoidal deformations (•), com- 
pared with the values from the mechanical routed (□), and 
the data of Thomson et aZ.— (A), El Bardouni et al^ (0), 
and Schrader et al^ (continuous curve); the planar limit is 
shown dotted. The predictions of FMT— are depicted in the 
inset (dashed curve). 
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tensor route (e.g., the recent MD data of Vrabec et al*&). 
This is in line with the concerns of Schofield and Hender- 
son^, and others ^ 20 ' 29 about the inadequacy of the me- 
chanical route for very small systems. For drops larger 
than R e ~ 8 we observe values of the surface tension 
which appear to be slightly larger than the planar limit, 
7(i?) > 7oo; because the tension has to converge to 
when R e — > 00 this suggests a non-monotonic behaviour 
of the tension with increasing curvature, and a corre- 
sponding weak maximum. Our values are consistent with 
the data point reported by El Bardouni et a/.— estimated 
from the surface free energy change, and with the small 
maximum observed by Schrader et al. 22 using a Landau 
free energy approach in the canonical ensemble (though 
the authors do not comment explicitly on this point). 
The calculations of the tension of curved interfaces from 
curvature corrections, SGT and DFT (which have been 
brought into question because of their failure to repro- 
duce existing simulation data) are now supported by our 
data. In the inset of Fig. 3 we compare the TA data 
for the surface tension of drops with those from a non- 
local DFT using fundamental measure theory (FMT)ii; 
a maximum is predicted with FMT at R e ~ 10. 

Three main conclusions can be gleaned from our study. 
Firstly, there is clearly a large fluctuation contribution 
to the interfacial tension of nanoscale spherical drops 
(and most likely other curved surfaces) in addition to the 
underlying first-order (mechanical) contribution which 
fully describes the planar interface. Such contributions 
from fluctuations in the energy are not found in the pla- 
nar limit to any significant degree. As a consequence, 
our results do not support the validity of a mechani- 



cal (pressure-tensor) route to the interfacial tension for 
surfaces of high curvature such as small drops. This 
is in line with the warning of Blokhuis and Bedeaux— 
that the use of a mechanical approach in this context 
"is still a matter of concern" and that it is "advisable 
not to use the pressure-tensor whenever this can be 
avoided" . Our data are not consistent with the mono- 
tonic dependence of the surface tension with curvature 
obtained from a mechanical treatment. As well as con- 
tributions in (x a (dU/dx a )), the correct "virial" expres- 
sion for the surface tension would have to contain terms 
in averages of the type (x a (dU /dx a )){xp {dU/dxp)) and 
(x a xp (dU/dx a ) (dU/dxp)), which would involve up to 
four-body correlations for pair-wise additive potentials. 
This suggests that there are additional contributions to 
the change in the entropy due to the deformation of small 
drops involving quadratic terms in ALU (At/ 2 ), (AU) 2 , 
(UAU 2 ), (AU)(UAU), (U){AU 2 ), and {U){AU) 2 . As 
a final point, the rise of the surface tension above that of 
the planar limit after a certain drop size would be con- 
sistent with a negative Tolman length. Our data for the 
larger drops suggests a value of<$oo ~ — 0.2 ±0.3. Though 
the statistical uncertainty is large, our finding supports 
the exact mean-field predictions for the penetrable sphere 
models, and is consistent with the latest accurate value 
of —0.10 ± 0.02 determined from the Laplace relation for 
a large N — 100, 000 particle system^. 
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